/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name : f_exp.c
 *  Function  :
 *        f_exp function returns an exponential value powed
 *        by given value.
 *  Revision  :
 *      2003/04/11    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/
#include <f_math.h>
#include <f_smcvals.h>


// ALGORITHM
//
//  exponential function can be showed by continued fraction expantion as below:
//  
//
//                             2x
//  exp( x ) = 1 + ------------------------------------         (a)
//                                 x^2
//                   2 - x + --------------------------
//                                      x^2
//                              6 + -------------------
//                                           x^2
//                                    10 + ------------
//                                          14 + ....
//
//  If you don't know it, search about "continued fraction" on the internet. 
//  Notify that Many common books about mathmatics doesn't have its information.
//
//  Its convergence speed becomes higher under | x | < (1/2) * log2. We assume
//  like this:
//
//     x = t + k * log2                       (b)
//
//  Here, value k can have only integer value. value t has only fraction part, doesn't
//  have integer part, and its allowed range of value is t < (1/2) * log2. Then exp(x)
//  can be showed like below:
//  
//     exp(x) = exp( t ) * 2^k                (c)
//
//   Finially, we try to get t and k, and get exp(x)'s value by f_frexp function of 
//  ANSI function.
//   In order to get k's value from x, devive (b) by log2.
//
//     x/log2 = t/log2 + k                    (d)
//
//   Here, t has only fraction part, so integer part of x/log2 is equal to k. t gets from
//  below:
//
//      t = x - k * log2                      (e)
//
//   we can get exp(t) of (c) by using (a), then get exp(x)'s value.
//
//    about more information, See below:
//      "Software technology13: most new algorithm dictionary by C"
//           written in Japanese by Haruhiko Okumura, 
//           published by gijutsu hyouron sya, 1991.5.10 third published
//
//
//   I added more flavors:
//       when x = 0,               exp( x ) = 1.0
//       when x < 0,               exp( x ) = 1/exp(-x)
//       when x >=87.3365484518 ,  exp( x ) = INF
//       when x < 0.01, try to multiply x = x* 2.0 , 5 times at most
//                                    ( only using double )
//
//  The larger x becomes, the worse f_exp(x) becomes. If you want to avoid it,
//  define F_MATH_DBL macro and re-compile f_math.lib. Then f_exp uses double 
//  calculation at some point.
//
//
//  the architecture of float floating point
//
//   0 30         23 22                0  32 bit
//   -----------------------------------
//  | |   exponent  |  fraction         |
//   -----------------------------------
//
//  |               |                   |
//  |    8 bits     |    23 bits        |
//
//         bit    31         sign bit        (  1 bit  )
//             30 - 23      exponent part    (  8 bits )
//             22 -  0      fraction part    ( 23 bits )
//
//




//#define F_MATH_DBL
#define LOG2 0.69314718055994530942f
#define N 22

float f_exp( float x )
{
	int i,k;
	float x2,w;
	unsigned long ulNaNData;

	// double using....
#ifdef F_MATH_DBL
	double sdX;
#endif

	int l;
	int errno_buf;
	FLT_LNG sfTemp1;


	sfTemp1._F = x;

	// exp( 0.0 ) = 1.0
	if ( (sfTemp1.st._LL & 0x7fffffff ) == 0 ) {
		sfTemp1._F = 1.0f;
		return sfTemp1._F;
	}

	if ( (sfTemp1.st._LL & 0x80000000) != 0x0 ) {

		errno_buf = errno;
		sfTemp1.st._LL &= 0x7fffffff;
		sfTemp1._F = f_exp( sfTemp1._F );

		if ( sfTemp1.st._LL == f_P_INF.st._LL ) {
			sfTemp1._F= 0.0f;
			errno = errno_buf;
		} else {

			/* If NaN, returns the returned raw value */
			ulNaNData = f_NAN.st._LL;
			if ( ( sfTemp1.st._LL & ulNaNData ) != ulNaNData ) {
				sfTemp1._F = 1 / sfTemp1._F;
			}

		}
		return sfTemp1._F;
	}

	/* below number is -87.3365484518.... */
	/* exp( -87. ... ) = float minimum number */
	if( sfTemp1.st._LL >= 0x42aeac50 ){

		if (  ( sfTemp1.st._LL & 0x7fc00000 ) == f_NAN.st._LL ) {
			/* Is NAN ? */
			errno = EDOM;
			return f_NAN._F;
		}
		errno = ERANGE;				// value check
		sfTemp1._F = f_P_INF._F;
		return	sfTemp1._F;

	}

	//  fraction part are added as 1.
	k = (int)(x / LOG2 + 0.5 );

	// double using....
#ifdef F_MATH_DBL
	sdX = (double)x;
	sdX -= (double)k * (double)LOG2;
	x = (float)sdX;
#else
	// float using....
	x -= k * LOG2;
#endif

	l = 0;

	//// double using...
#ifdef F_MATH_DBL
	sfTemp1._F = x;
	for ( l = 0 ; sfTemp1.st._LL <= 0x3a83126f, l < 6 ; l++ ) {
		sfTemp1._F = sfTemp1._F * 2.0f;
	}
	x = sfTemp1._F;
#endif

	x2 = x * x; w = x2 / N ;
	for( i = N -4; i >= 6 ; i -= 4 ) w = x2 / ( w + i );

	return f_ldexp( ( 2 + w + x ) / ( 2 + w - x), k +l );
}
